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Preliminary  Re marks 

This  report  can  be  said  to  have  two  parts.  The  first  part  is  an 
elucidation  of  some  notions  and  principles  assembled  by  the  senior  writer  to 
expedite  the  calculation  of  the  lag  function. 

The  second  part  consists  of  detailed  instructions  on  how  to  program  the 
concepts  presented  in  the  first  part.  It  was  written  by  the  senior  writer 
with  a  view  to  expediting  the  understanding  of  an  unknown  program  writer  of 
unknown  capacities.  As  a  consequence,  it  is  perhaps  too  elementary  in 
character  to  appear  here.  However,  an  instructional  precept  of  some  merit  is 
that  it  is  better  to  over-  rather  than  under-simpli f v .  Time  loss  to  the 
knowledgeable  should  be  negligible  if  scanning  is  used  judiciously. 


at  end  (if  paper. 
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The  j  in  the  exponential  function  is  the  imaginary  unit  and,  since 

x+iy  x  x 

e  =  e  cos  y  +  i  e  sm  y  , 

one  has 

2jkfq2-t2]  ,  2  2  s  .  .  ,  2  2  s 

e  ;  =  cos  2k  (_q  -  t  J  +  l  sin  2k (q  -  t  J 

=  cos  2kq2  cos  2kt2  +  sin  2kq2  sin  2kt2 

2  2  2  2 
+  i[sin  2kq  cos  2kt  -  cos  2kq  sin  2kt  ] 

Complex  arithmetric  is  ruled  out  of  the  calculation  in  order  to  preserve 

2  2  s 

desirable  attributes  of  the  expanded  form  of  e'2^^  1  so  that  the  two 

integrals  become  eight  integrals  with  each  integrand  containing  a 
trigonometric  product  as  a  factor. 

Let  the  q-axis  be  vertical  and  the  t-axis  horizontal.  Then  each  of  the 
products  will  be  seen  to  vanish  on  a  set  of  horizontal  lines  and  a  set  of 
vertical  lines  that  subdivide  the  two  areas  of  integration,  infinite  strips 
bounded  by  q  =  0  and  q  =  1,  into  rectangles  in  which  the  trigonometric  factor, 
moving  either  horizontally  or  vertically,  alternates  in  sign.  As  a 
consequence,  if  the  integrals  are  considered  as  composed  of  the  sum  of  the 
integrals  on  the  rectangles,  then  each  double  integral  is  equivalent  to  the 
sum  of  the  sums  of  a  finite  set  of  alternating  series  infinite  in  the  t- 
direction.  It  follows  that,  if  each  integral  is  evaluated  by  the  rectangles 
in  a  row,  row  by  row,  the  question  of  "far  enough"  becomes  the  question  of 
far  enough  for  the  infinite  series.  Since  the  algebraic  component  in  the 
integrand  of  each  integral  does  not  change  sign  in  the  area  of  integration, 
the  series  are  strictly  alternating  and  so  the  far-enough  question  is  easily 


answered . 
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The  problem  of  evaluating  the  lag  function  thus  reduces  to  evaluating  a 
number  of  double  integrals,  each  with  an  integrand  well  behaved  in  the 
rectangle  of  definition.  Evaluating  a  multiple  integral  by  iteration,  6  has 
a  cyclical  character.  Algorithms  with  cyclical  properties  tend,  in  program 
form,  to  be  relatively  simple,  which  suggests  that  iteration  is  the  method  of 
choice.  To  write  a  program  in  iterated  form  for  a  multiple  integral,  one 
first  chooses  a  one-dimensional  quadrature.  This  quadrature  's  then  expanded 
to  two  (or  more)  dimensions  by  writing  its  quantities,  except  for  an  obvious 
few,  with  one  additional  independent  variable,  the  additional  variable 
indicating  a  level  of  integration  in  the  Iteration. 

Controlling  Error  in  Numerical  Integration 

As  anyone  at  all  experienced  in  the  art  of  computing  well  knows,  a  prime 
difficulty  in  the  computing  of  tables  is  the  problem  of  keeping  the  incidence 
of  unacceptable  errors  sufficiently  small,  for  errors  seem  to  creep  in  with  m 
inevitability  that  sometimes  seems  to  suggest  the  guidance  of  a  malign  hand. 
The  true  cause,  of  course,  is  the  many  points  of  possible  entry,  ranging  from 
faulty  computing  procedures  to  a  simple  inversion  of  digits  in  those  inputs 
and  outputs  handled  by  a  human  hand.  The  best  of  all  checks,  one  seldom 
available,  is  the  existence  of  a  relationship  between  the  final  outputs  which 
is  known  to  hold.  Unfortunately,  no  such  check  seems  to  exist  for  the  final 
outputs  of  the  lag  function  except  the  one  always  possible  for  results  in 
tabular  form,  namely,  finite  differencing.  The  values  given  for  the  lag 
function  were  subjected  to  this  test,  although  the  test's  efficacy  was 
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diminished  some  by  the  necessity  o£  using  mostly  divided  differences  rather 
than  ordinary.  Divided  differences  tend  not  to  indicate  a  faulty  entry  as 
precisely  as  ordinary  differences  do. 

To  make  up  for  the  Lack  of  a  dependable  final  check,  considerable  care 
was  taken  to  ensure  that  the  accuracy  of  each  of  the  individual  numerical 
integrations  was  good  to  three,  if  not  four,  significant  figures.  In  general, 
there  is  no  means  of  estimating  an  error  in  an  evaluation  by  a  particular 
quadrature  formula  other  than  comparing  it  with  another  evaluation  using 
another  quadrature  differing  in  some  significant  respect.  Some  comparisons, 
however,  do  not  always  give  reliable  estimates  because  of  a  cancellation  of 
errors . 

A  favored  procedure  for  numerical  evaluation  of  a  definite  integral 
begins  by  dividing  the  interval  of  integration  into  subintervals.  The 
subdividing  is  likely  to  be  guided  to  some  degree  by  the  perceived  character 
of  the  integrand  and  by  the  choice  of  the  quadrature  formula  which  is  to  be 
used  in  the  subintervals.  The  chosen  quadrature  obtains  from  the  set  of 
subintervals,  a  set  of  subintegrals  which  when  summed  is  taken  as  a  first 
approximation.  A  second  approximation  is  then  formed  by  using  the  same 
quadrature  on  the  set  of  subintervals  that  come  from  halving  the  lengths  of 
the  subintervals  in  the  first  set.  The  two  evaluations  are  compared  by 
taking  their  difference.  If  the  difference  is  sufficiently  small,  the 
second  approximation  is  then  accepted  as  final. 

Suppose  the  quadrature  chosen  is  one  of  the  Newton-Cotes  kind,  say, 
Simpson's  rule.  Consider  the  subapproxi.nat  ion  obtained  on  one  subinterval  of 
the  first  set.  It  can  be  written  in  the  form 
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where  h  is  the  interval  length  and  yq  and  V4  are  ordinates  at  its  end  points 
and  y2  Ht  its  midpoint.  The  interval  length  for  the  matching  two  subintervals 
of  the  second  set  is  h/2  and  their  quadrature  formula  becomes 

XT  iy 0  +  4yt  +  2y2  +  4y2  +  y4)  • 

The  difference  between  the  two  subresults  will  be  seen  as  h  multiplied  into 
a  fourth  order  finite  difference:  hA4y/12. 

That  the  difference  between  the  two  quadratures  is  found  to  be  a  constant 
multiple  of  a  fourth  order  finite  difference  is  not  entirely  happenstance. 

For,  as  will  be  shown,  the  difference  between  any  two  quadratures  on  the  same 
interval  is  always  a  linear  combination  of  finite  differences  of  the  order  of 
the  lower  of  the  two  quadratures.  In  the  case  just  considered,  both 
quadratures  are  of  the  fourth  order,  and  accordingly,  their  difference  must  be 
a  linear  combination  of  fourth  order  finite  differences,  which  it  is,  though 
it  consists  of  only  one  term. 

A  more  interesting  example  is  had  in  the  difference  between  Simpson's 
rule  and  the  four  term  Newton-Cotes  formula: 

h(y0  +  'Jy1  +  3y2  +  y3)/S  . 

To  have  the  seven  abscissas  of  the  two  formulas  fall  evenly  in  a  common 
interval,  it  is  necessary  to  subdivide  the  interval  (0,h)  into  six  equal 
parts.  The  seven  abscissas  are  then  given  by  the  relation 

x^  =  kh/b,  k  =  0,1, 2,... 6 


It  follows  that  Simpson's  rule  takes  the  form 
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h (y0  +  4y3  +  y6}/&  . 

the  four  term  formula,  the  form 

h(y0  +  3y2  +  3y4  f  y.'J,,!3 

and  the  difference  may  be  written  as 

ii(y6  -  9y4  +  16y3  -  9y>2  +  yQ)/43  . 

The  parenthetical  quantity  with  y3  and  y^,  added  and  subtracted,  is 
equivalent  to 

(y 6  ~  y$J  +  (y5  -  y4 )  -  3 (y4  -  y3)  f  s(y3  -  y2)  *  (va  '  yi^  f  -yi  "  yiv 

or,  in  terms  of  first  order  finite  differences,  to 

^5  +  Ay^  -  8Ay3  +  3  Ay.,  -  Ayt  -  Ay()  . 

This  conversion  of  a  linear  combination  of  finite  differences  of  order  zero 
to  one  of  the  first  order  is  accomplished  by  grouping.  Additional  groupings, 
three  in  number,  reduce  the  initial  linear  combination  to  one  of  the  fourth 
order*:  A4yo  +  4A4y^  +  A4yQ  , 

•Jf 

The  sane  result  can  be  obtained  more  easily  by  a  parallel  to  synthetic 
division.  Consider  the  following  array  which  is  like  that  of  synthetic 
di vi s ion . 

I  0  -9  lb  -9  0  1  U 

1  1-8  3-1-1 

11-8  8  -1  -1  |_1_ 

1  2  -6  2  1 

12- 6  2  1  |J_ 

l  3-3-1 

13- 3-1  |J_ 

1  4  1 


1  4 
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The  first  course  of  numbers  in  this  display  are  the  numerical 
coefficients  in  the  linear  combination  of  finite  differences  of  order  zero; 
the  second,  the  coefficients  in  the  linear  combination  of  first  order 
differences  and  so  on,  course  by  course.  The  last  fully  interpreted  is 

4  4  4 

A  V 2  +  -*  A  y  j  +  A  y0  . 

a  linear  combination  of  the  kind  designated  by  the  principle  stated  above. 
The  full  expression  for  the  difference  between  the  two  original  quadratures 
becomes 

'4  4  4  \ 

hi  A  w  y  +  A  y^  +  A  y  j/48 

The  principle  gives  only  the  order  of  the  finite  differences  in  the 
final  linear  combinat ions .  It  does  not  give  the  coefficients  or  the  total 
number  of  the  differences.  The  latter  quantity,  however,  can  be  deduced 
from  the  total  number  of  abscissas  in  the  difference  of  the  quadratures. 

The  number  of  abscissas  involved  in  the  case  using  Simpson's  rule  is  five. 
On  five  points  only  one  fourth  order  finite  difference  is  possible.  In 
the  second  case  with  the  addition  of  abscissas  X9  and  X4 ,  the  total  is 
seven.  Without  going  through  the  whole  procedure,  one  knows  that  only 
three  finite  fourth  order  differences  are  possible.  So  all  that  is 
Learned  from  the  full  procedure  is  that  the  multipliers  (1,4,1),  each 
divided  by  48,  will  appear  in  the  final  result. 

Both  error  estimates  in  the  two  illustrative  examples  are  combinations 
of  ordinary  fourth  order  finite  differences.  So,  if  the  integrands  in  the 
subintervals  have  small  variation,  the  individual  errors,  anil  also  the  sum 


\ 


-11- 


27  March  1981 
GHP'.SJG: Ihz 


sum  of  their  absolute  values  over  the  subintervals ,  can  be  expected  to 
diminish  with  some  rapidity  as  the  subintervals  are  halved.  In  the 
instance  of  the  second  example,  however,  cancellation  between  the  three 
finite  differences  may  indicate  an  error  appreciably  less  than  actual  and 
allow  acceptance  of  an  incorrect  result.  This  nay  seem  very  improbable 
but  the  writer's  experience  suggests  that,  when  it  comes  to  errors, 
Murphy's  law  always  should  be  remembered. 

The  Newton-Cotes  formulas  with  their  even  spacing  often  have  a 
desirable  simplicity  in  some  applications  of  numerical  integration.  They 
are,  as  Is  well  known,  much  less  powerful  than  those  of  the  Gaussian 
kind.  The  calculation  of  the  lag  function,  though  not  awe  inspiring 
by  today’s  standards,  is  still  large  enough  to  warrant  using  Gaussian 
quadratures . 

The  principle,  whose  validity  was  established  in  two  instances  with 
quadratures  of  the  Newton-Cotes  kind  also  holds  for  the  Gaussian,  the 
Radau,  the  Lobatto  quadratures  and  others,  but  because  of  the  uneven 
spacing  of  their  abscissas,  ordinary  finite  differences  must  be  replaced 
by  divided.  Divided  differences  do  not  yield  so  readily  to  manipulation. 
A  proof  of  general  validity  must  be  found. 

One  illustrative  example  tells  in  essence  the  whole  story.  Suppose 
the  two  quadratures  are  the  two  term  and  the  three  term  Gaussian 
quadratures  on  the  same  interval,  bet  the  two  quadratures  be  represented 


by  the  quantities 


-12-  27  March  1981 

GHP-.S.JG:  lhz 

VO  +  c2y>  +  C4y4  • 

where  Che  c's  are  weights  and  the  v's  ordinates.  Consider  their  difference 

c0-v0  “  Vl  +  c2-v2  ~  °3'  '3  +  c4y4  • 

The  two  term  Gaussian  gives  the  exact  value  of  the  integral  of  any  polynomial 
of  third  degree,  the  three  term  Gaussian,  of  any  polynomial  of  the  fifth 
degree.  Let  the  polynomials  be  the  first  four  powers  of  x,  namely,  1,  x, 
x-  and  xJ.  Both  formulas  give  the  exact  value  of  the  integrals  for  these 
powers  of  x.  When  each  power  is  used  in  succession  to  supply  values  for  the 
ordinates  in  the  difference,  one  obtains  a  system  of  four  homogeneous 
equations  in  the  c's: 

4 

~  c.xk  =  0  ,  k  =  0,1 . 4  . 

J  J 

.1=0 

The  rank  of  the  matrix  of  the  system  is  three.  Therefore,  a  solution  of 
the  four-equation  system  exists,  other  than  c^  =  Cj  =  e.,  =  c^  =  c^.  Any 
constant  multiple  of  this  solution  is  also  a  solution.  The  fourth  order 
finite  difference  on  the  five  points  fx^y^),  (x^vj,  l x  -j » >'  3 ) 

and  fx, ,v.  1  is  also  a  linear  combination  of  ordinates  and  also  must  vanish 
for  the  four  powers  of  x.  In  view  of  the  generality  of  this  last  proof  of 
the  principle  stated  above,  it  must  be  true  for  any  two  quadrature 


formulas 
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A S access  ion  of  Exp  Lan a t o rv  Alport t  h m s P r e 1 i m i n ary 

to  the  Working  Algorithms 


Although  the  senior  author,  over  the  years,  has  formulated  a  number  of 
algorithms  for  expression  and  calculation  in  the  language  of  FORTRAN,  he  lias 
never  written  a  program  in  the  language,  only  through  the  services  of  an 
intermediary.  Illiteracy  in  FORTRAN  and,  for  that  matter,  all  computer 
languages  persisted  until  the  appearance  of  JOSS'',  a  language  so  obvious 
and  consistent  in  terminology  and  forms  that  to  Learn  it  is  a  negligible 
task  and  to  write  a  program  in  it,  a  pleasure. 

In  the  case  of  the  lag  function,  several  algorithms  for  approximating 
its  definite,  integrals  were  written  in  JOSS  and  checked  for  faults.  All 
use  basic  Gaussian  quadratures:  one  of  two  terms,  one  of  three,  and  one 
of  four.  The  K-th  abscissa  of  the  quadrature  of  J  terms  is  represented 
by  D(J,K)  and  the  K-th  weight  by  C(P,K).  For  the  interval  (0,1),  the  values 
of  the  C's  are: 


C( 2  ,  l )  =  .5 

C( 3  ,  L )  =  .2777  ... 

C( 3 , 3)  =  C( 3 , 1 ) 
C(4 ,2)  =  .326072577431273 
C( 4 , 3)  =  C(4 ,2) 


C( 2 , 2 )  =  C( 2,1) 

C(  3 ,2)  =  .444  ... 

CC  4,1)  =  .173927422568727 
C(  4 , 4  )  =  C( 4 , 1 ) 


and  the  values  of  the  D’s  are: 

0(2,1)  =  .5  -V3/b 
0(3,1)  =  .112701665379258 

D( 3 , 3)  =  1  -  0(3,1) 
D ( 4 ,  l )  =  .069431344202973 

D(4,3)  =  1  -  D(4 ,2) 


D( 2 , 2 )  =  1  -  0(2, l) 

0(3,2)  =  .5 

0(4,2)  =  .330009478207572 
0(4,4)  =  1  -  (4,1) 


An  acronym  for  "Johnniac  Open  Shop  System" , 
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The  first  ot  the  explanatory  algorithms,  explanatory  with  regard  to  JOSS 
and  to  the  working  algorithms,  follows: 

Algorithm  1 


1.12 

Set  11  =  8  -  A 

1  . 18 

Set  J  =  1 

i  .3 

Set  V(J)  =  0 

1.38 

Set  R  =  0 

1 .4 

Set  J  =  J  +  1 

1 .42 

Set  V(J)  =  0 

1.5 

Set  K.  =  K  +  1 

1.52 

Set  X  =  A  +  D( J 

,K)*H 

1.56 

Set  V(J)  =  V(J) 

+  V  (  J  , 

,R)*H*F(X) 

1 .58 

To  Step  L . 48  if 

K  <  J 

l  .62 

To  Step  L.94  if 

t  V  ( J  ) 

-  V  (  J  - 1 )  1 

1 .64 

To  Step  1.38  if 

J  <  4 

1  .94 

Type  V( J ) . 

This  program,  which  is  capable  of  approximating  the  definite  integral 
of  F(X)  on  the  interval  (A,B),  exhibits  several  essential  characteristics 
of  the  JOSS  language:  Programs  are  written  in  parts,  where  a  part  is  similar 
to  a  FORTRAN  subroutine  with  a  return.  A  part  is  made  up  of  a  succession  of 

numbered  commands.  The  number  on  a  command  is  the  sum  of  an  integer,  which 

does  not  differ  for  commands  from  the  same  part  and  a  decimal  fraction. 

The  commands  are  grammatical  to  the  extent  that  the  first  letter  of  the  first 
word,  always  a  verb  in  the  imperative  mood"  ,  is  capitalized  and  the  end  of 
the  command  is  signaled  by  a  period. 

The  flow  of  the  program  is  easily  grasped.  The  body  of  the  program 
implements  the  three  quadratures  one  by  one  in  a  evclical  fashion.  The 

quantitv  K  in  Step  1.62  specifies  an  allowable  error  and ,  barring  the 


*See  Appendix  A  for  FORTRAN  equivalent. 
**. 


The  verb  "go"  in  Steps  1.68,  1.62  and  1.64. 
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aruisa.il,  if  V(J)  is  accepted  for  J  =  2,  3  or  4  ,  the  error  in  the  approximation 
obtained  is  presumably  less  than  or  equal  to  K.  If  V(2)  is  accepted,  it 
follows,  because  Step  1.3  sets  V ( 1 )  to  zero,  that  the  approximation  is  less 
than  K. 

When  V(4)  is  not  accepted  by  Step  l.b2,  the  control  is  not  sent  to 

Step  1.38,  but  passes  on  to  Step  1.94  for  a  print  out.  Inasmuch  as  the 

quadrature  of  four  terms  yieLds  the  correct  result  for  a  polynomial 

integrand  of  the  seventh  degree,  the  small  program  may  well  be  adequate 

for  the  integrands  it  will  be  called  on  to  evaluate.  After  all,  the 

integrands  it  needs  to  approximate  do  not  necessariLv  behave  in  an  outrageous 

fashion.  To  cover  the  contingency  that  some  will  be  difficult,  the  small 

program  is  so  modified  in  Algorithm  2  (which  follows)  that  if  V(4)  is  not 

accepted,  the  interval  (A,B)  is  halved  and  each  half  attacked  separately. 

If  the  interval  of  integration  is  to  be  halved  and  rehalved,  there  are 

many  ways  to  go  about  it.  One  way  is  presented  in  the  following  algorithm: 

Algorithm  2'* 

l.l  Set  0  =  l 
1.12  Set  11(0)  =  B  -  A 
1.14  Set  U(0)  =  A 
1.1b  Set  S  =  0 
1.18  Set  J  =  1 

1.3  Set  V(J)  =  0 
1.38  Set  J  =  J  +  1 

1.4  Set  K  =  0 
1.42  Set  V(J)  =  0 

1.5  Set  rl  =  K.  +  1 

1.52  Set  X  =  i;(U)  +  D(J,K)*H(U) 

1.5b  Set  V(J)  =  V(J)  +  C(.J,K)*M(0)*F(X) 

1.58  To  Step  1.5  if  K  <  J 

1  .82  To  Step  1.18  if  |V(.J)  -  V(J-l)  <  K 

1.64  To  Step  1.38  if  .1  <  4 


*See  Appendix  A  for  FORTRAN  equivalent. 
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1.6b  11(0)  =  H  ( 0 )  /  2 

1.68  Set  0  =0+1 

1.7  Set  H(0)  =  11(0-1 ) 

1.72  Set  U(0)  =  U(0-1)  +  11(0) 

1.76  To  Step  1.18 
1.84  Set  S  =  S  +  V(J) 

1.86  To  Step  1.94  if  0  =  1 
1.88  Set  0  =  0  -  l 
1.9  To  Step  1.18 
1.94  Type  S. 

In  subdividing  an  interval,  the  symboL  0  is  introduced  whose  largest 
value  at  any  time  is  equal  to  the  number  of  intervals  available  to  the 
algorithm  for  processing.  All  intervals  are  tagged  with  an  0-value.  At 
the  beginning  there  is  only  one  interval,  [4,B].  So  0  =  1.  If  the  three 
quadratures  are  unable  to  reach  an  acceptable  result  on  this  interval,  it 
is  divided  into  two  equal  parts,  the  part  nearest  B  is  given  the  value 
0=1,  the  other  part  retains  the  value  1.  The  attack  shifts  to  the 
interval  with  0=1.  If  the  attack  is  successful,  the  attack  moves  to 
0=1.  If  unsuccessful,  0  =  2  is  divided  and  the  two  parts  tagged  with 
the  cardinal  numbers  0=2  and  0=3.  Always  the  interval  with  the 
highest  0-value  is  the  one  in  process  until  ultimately  success  comes  to 
the  interval  0=1. 

Algorithm  2,  it  will  be  noticed,  is  flawed  with  respect  to  the  control 
of  error.  Subdivision  of  intervals  leads  to  an  accumulation  of  error  that 
may  exceed  the  value  of  F  specified  in  Step  1.62.  Although  t lac  algorithm 
can  be  modified  so  as  to  insure  an  error  less  than  F,  the  simplest  and  a 
likely  adequate  modification  is  a  carefuLly  chosen  decrease  in  the  value 
of  E,  particularly  In  view  of  the  fact  that  there  are  a  number  of  other 
sources  of  error  to  protect  against  on  the  way  to  the  value  of  B(k). 

The  second  algorithm  wiil  now  be  used  to  illustrate  the  conversion 
of  an  integration  in  one  dimension  to  an  iterated  integral  in  two 


dimensions . 
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A_l;V>J^i  thm  _3'c 

1.1  Set  0(1)  =  1 

1.12  Set  Hi  1,0(1)]  =  B(l)  -  A( I ) 

1.14  Set  U [  1,0(1)]  =  A( 1 ) 

1.16  Set  S(l)  =  0 

1.13  Set  J(I)  =  1 

1.3  Set  VU,.J(I)]  =  U 

1.38  Set  JQ)  =  J(I)  +  1 

1.4  Set  K(i)  =  0 

1.42  Set  V [ I , J ( I ) ]  =  0 

1.5  Set  K.(  I)  =  K(I)  +  1 

1.52  Set  X(I)  =  U[ 1,0(1)]  +  D [ J ( 1 ) ,  K(  I )  ]  *  H  [ 1,0(1) ] 

1.54  Do  Part  2  for  I  =  2  if  I  =  l 

1.56  Set  V[I,J(I)]  =  V[ I, J( I) ]  +  C[J(L),K(I)]*ll[I,0(I)]*G(I) 

1.58  To  Step  1.5  if  K(I)  <  J(I) 

1.62  To  Step  1.84  if  |  V  [  I ,  J  (  I )  ]  -  V[I,.J(L)  -  1]|  <  K(L) 

1 .64  To  Step  1.38  if  J(I)  <4 

1.66  Set  Hi  1 ,0( I) ]  =  H[I,0(I)]/2 
1.68  Set  0(1)  =  0(1)  +  1 

1.7  Set  H [  1 , 0 (  I )  ]  =  H]  1 ,0(  I)  -  1] 

1.72  Set  U[ I ,0( I) ]  =  Ull, 0(1)  -  1]  +  HlI,0(I)] 

1.76  To  Step  1.18 

1.84  Set  S(I)  =  S(l)  +  V [  I ,  J(  I ) ] 

1.86  To  Step  1.94  if  0(1)  =  1 
1.88  Set  0(1)  =  0(1)  -  1 
1.9  To  Step  1.18 
1.94  Type  S(l)  if  I  =  1 

2.3  Do  Part  1 

2.7  Set  l  =  1 

Except  for  the  index  I,  the  last  algorithm  and  the  preceding  one  are 
identical  for  the  first  11  steps.  Step  1.54  has  no  counterpart.  The 
counterpart,  if  it  existed,  would  fall  between  Steps  11  and  12.  Step  12 
and  those  following  closely  parallel  Step  1.56  and  following.  The  function 
G(I)  of  Step  1.56  is  defined  bv  the  relations 

I  S(2)  ,  if  1  =  1  , 

G(I)  =  1 

|  F l X( 1) , X( 2) ]  ,  if  1  =  2  . 


*See  Appendix  A  for  FORTRAN  equivalent. 
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The  vaLue  of  G(L)  for  L  =  2,  identities  the  iterated  interal  as  c 

8(1)  8(2) 

J  dx  /  F(x,v)  dy 

A(l)  A(2) 

it  helps  in  understanding  the  two-dimensional  algorithm  to  call  to  mind 
the  principles  of  the  iterated  integral.  As  an  example,  take  the  iterated 
integral  displayed  above.  The  last  integration  is  a  definite  integral  with 
respect  to  x  having  as  its  integrand  the  result  of  a  definite  integral  with 
respect  to  y  the  latter  having  F(x,y)  as  an  integrand  with  x  treated  as  a 
constant.  If  the  steps  in  the  algorithm  are  examined  careful Ly,  it  is  seen 
that  the  algorithm  in  its  own  way  follows  the  same  path.  The  algorithm  starts 
with  a  quadrature  with  respect  to  x,  or  rather  X( 1 ) .  The  first  abscissa  is 
the  first  abscissa  of  the  two-point  Gaussian  quadrature.  But  G(I)  in 
Step  1.56  requires  S(2),  which  is  the  resuLt  of  a  quadrature  with  respect  to 
y  [or  X ( 2 ) 1  with  the  integrand  F[ X( 1 ) , X( 2 ) ] ,  the  variable  X(  1 )  being  held 
constant.  Hence,  the  quadrature  with  respect  to  X( 1 )  must  be  temporarily 
abandoned.  This  is  accomplished  by  Step  1.54  which  sends  the  control  to 
Part  2  with  1=2.  The  first  step  of  Part  2,  Step  2.3,  sends  the  control 
back  to  Part  1,  i.e.,  Step  1.1.  The  quadrature  with  respect  to  X ( 2 )  then 
begins  and  continues  without  interruption  until  completed.  Control  then 
returns  to  Step  2.7,  which  sets  1=1,  completing  Part  2.  Control  returns 
to  Step  1.56  with  the  required  values  of  S(  2 )  now  available.  Hence  the 
quadrature  with  respect  to  X( 1 )  can  continue  until  shortly  Step  1  .  34  is 
encountered  again,  and  so  on. 

K( 1 )  and  2(2)  should  not  he  given  the  same  value  because  the  integrand 
for  the  X( 1 )-quudrature  has  a  component  from  the  X( 1 )-quad r  itur"  and  so 
many  have  a  scatter  is  Large  as  K(2).  Finite  differences  ire  well  known 
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for  their  sensitivity  to  scatter.  Consequently,  the  programmed  check  for  the 
X( l )-quadrature  tends  to  magnify  the  apparent  estimate,  making  acceptance 
difficult,  or  even  impossible.  Experience  indicates  that  E(l)  and  E(2)  should 
satisfy  the  relation, 

E( 1 )  <  E( 2 ) / 10  . 


The  Algorithm  for  the  Integral 


/  sin 
0 


f  sin  2kq 


2  ^  fl~9~ 
q+t 


dq 


The  second  integral  in  the  definition  of  B(k),  unlike  the  first,  is  free 

of  such  nuisances  as  discontinuities  in  its  integrand  and,  if  one  selects 

2  2N 

sin  2kt”  sin  2kq^  from  the  four  components  of  e^2*^  1  as  the  trigonometric 

factor,  one  has  the  simplest  of  the  eight  that  are  implicit  in  B(k)'s 
definition.  The  quantity  sin  2kx~  vanishes  when  x  =  ~\J  nrt/(  2k)  and 
n  =  0,1,2,...  .  Hence  the  area  of  integration  is  divided  into  rectangular 
sub-areas  by  the  lines 


9 


t 


where  p  =  m/(2k),  n  =  0,1,2,...  and  n  =  0,1,2,...  .  The  integration  over  the 
infinite  strip  is  accomplished  in  a  piecewise  fashion.  Starting  with  the 
rectangle  common  to  the  first  row  and  first  column,  the  calculation  proceeds 
rectangle-by-rectangle  along  the  row  until  a  suitable  test  indicates  the 
remainder  for  the  infinite  alternating  series  is  sufficiently  small.  The 
next  row  is  then  given  the  same  treatment  and  finally  the  finite  series  of 
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2  •> 

sums  is  summed  to  give  the  value  of  the  sin  2kt“  sin  2kq~  component  of  the 
second  integral.  The  steps  in  the  algorithm,  their  explanation  delayed, 
follow. 

1 . 1  Set  p  =  tt /  (  2*k) 

1.16  Set  N( 2 )  =  1 

1.18  Set  A( 2)  =  0 

1.2  Set  NCI)  =  1 
1 .22  Set  A( 1)  =  0 
1.24  Set  T[NC2) ]  =  0 

1.26  Set  B( 2)  =  sqrt[N(2)*p] 

1.28  To  Step  1.32  if  B( 2 )  <  1 

1.3  B(  2)  =  1 

1.  32  Set  B(  1 )  =  sqrt  [  N'(  l  )*p] 

1.38  Oo  Part  2  for  I  =  1 

1.42  Set  T l N ( 2 ) 1  =  T [ N ( 2 ) ]  +  SC  1 ) 

1.44  To  Step  1.56  if  NCI)  <  i(l) 

1.46  To  Step  1.56  if  [  S ( 1 ) j  >  EC  3 ) 

1.48  To  Step  1.56  if  s*S(l)  >  0 

1.5  Set  Z(l)  =  ZC1)  +  1 

1.52  To  Step  1.66  if  Z(l)  >  iCD 

1.54  To  Step  1.58 

1.56  Set  7.(1)  =  0 

1 . 58  Set  s  =  S(  l ) 

1.60  Set  A( 1 )  =  B(l) 

1 .62  Set  N( l)  =  N( 1)  +  1 

1.64  To  Step  1.32 

1.66  To  Step  1.88  if  B(2)  =  1 

1.68  To  Step  1.8  if  N(2)  <  i( 2) 

1.7  To  Step  1.8  if  |  T [ N  C  2 ) ]  |  >  E(3) 

1.72  To  Step  1.8  if  T[N(2)]*T[N(2)  -  1]  >  0 

1.74  Set  7.(2)  =  Z(2)  +  1 

1.76  To  Step  1.88  if  Z(2)  >  i ( 2 ) 

1 .78  To  Step  1 .82 

1.8  Set  Z( 2)  =  0 

1.82  Set  A( 2)  =  8(2) 

1.84  Set  N( 2 )  =  N(2)  +  1 

1.86  To  Step  1.2 

1 .88  Print  output 

2.1  Set  0(1)  =  1 

2.12  Set  H[ i,0( I) ]  =  B( I)  -  A( I ) 

2.14  Set  Ul  1,0(01  =  A(l) 

2.16  Set  S(I)  =  0 

2.18  Set  J(I)  =  1 

2.3  Set  V( I,J(I)1  =  0 

2.38  Set  .1(1)  =  .1(0  +  1 
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2  .4 

Se  t  K  (  l  )  —  1 

2.42 

Set  V l I , -H  1  )  ]  = 

,i 

2 . 3 

Se  t  K(  t )  =  l\(  I  ) 

b  1 

2.32 

Set  K(i)  =  U  L  L  , 1 

KOI  +  U|J(0,K(l)l*!ll  I,  KOI 

2.34 

Do  Part  3  for  1 

=  2  if  I  =  l 

2.3p 

Set  V [ I , J ( l ) 1  = 

V|I,J(1)1  +  C|J(D,K(D!*il[l,u(I)l*d(D 

2.38 

To  Step  2.3  if 

K(  1)  x  J(  1) 

2 .  ti  2 

To  Step  2.84  if 

i v t i , . i ( o i  -  vo,j(o  -  m  s  s(o 

2.66 

Set  llll.fKOl  = 

ul  1,0(01/2 

2 .  b8 

Set  0(1)  =  o(i) 

+  l 

2.7 

Set  UU,0(1)]  = 

hi  1, 0(0  -  U 

2.72 

Set  Ul 1,0(1)]  - 

ill  ,0(0  -  1 1  +  111  L,0(  01 

2.76 

To  Step  2.18 

2.8 

Set  S(  [)  =  S(  1 ) 

+  v( I ,J( D  ! 

2.86 

Done  if  d( 1 )  - 

1 

2.88 

Set  0(0  =  "( O 

-  1 

2.9 

To  Step  2.18 

3.3 

Do  Part  2 

3.7 

i  =  l 

Inputs  tor  the  a!  torithm  ar  ■,  <>i  omirso ,  k,  tlio  weights  and  fnscissas 
for  the  throe  Gaussian  quadr  itures ,  aiul  the  oritioaL  quantities  K(  O  and 
H  ( — )  -  In  addition,  there  is  a  third  critical  quantity,  ';■.()),  and  two 

critical  integers  ill)  and  i(2).  i'heir  roles  are  made  clear  in  the 
foL Lowing  explanation  of  the  steps  of  the  algorithm. 

Part  2  carries  the  main  body  ot  the  eosnputat ion  and  can  be  said  ta  he 
served  by  Parts  1  and  3.  Steps  prior  to  Step  1 . IS  prepare  the  wav  for 
Part  2,  Step  1.38  initiates  the  computations  of  Part  2  and,  when  Part  2 
has  finished  with  a  rectangle,  steps  subsequent  to  Step  1 . 38  sum  the 
result,  select  the  next  rectangle,  if  another  is  required,  and  take  over 
tlie  functions  of  some  of  the  steps  prior  to  Step  1.38. 

The  quantities  N ( 1 )  and  N(23  in  Steps  l.lo  and  1.13  are  defined  bv 
stating  that  the  number  pair  | \(  1 )  ,  \‘(  2 )  ]  are  the  referene  numbers  lor 
the  rent  untie  common  to  the  Ml)^'  column  and  N(  2 ) c''  row. 

T[\(2)l  of  Step  1.2a  is  the  running  total,  iccumlated  rec t amt  1 e-hv- 


rectaintle  aloiut  the  N(2)c*1 


row . 
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Stops  1,44  through  l .  5h  express  properties  or  convergent  e.lternutins 
series.  It  will  be  recalled  that  when  a  series  readies  a  term,  after 
which  it  is  strictly  alternating  and  the  magnitude  of  each  ten  is  less 
than  the  procedi  hr  terra,  then  tire  magnitude  of  the  last  terra  calculated 
is  a  hound  on  the  remainder.  The  quantity  2(1)  in  Steps  1.3,  1.32  and 
l, 5b  is  a  count  of  the  number  of  times  in  succession  the  terms  iron  a 
row  have  met  these  conditions.  Step  1.44  does  not  allow  the  count  to 
start  until  N( 1 )  reaches  a  value  thought  to  be  large  enough  to  make 
acceptance  by  happenstance  unlikely.  The  character  of  the  integrand 
insures  that  such  a  value  exists.  Step  l.in  returns  the  count  to  zero, 
whenever  5(1),  the  magnitude  of  the  '.'(l)1-'1  terra,  is  not  loss  than  an 
appropriately  chosen  critical  quantity.  Step  1.48  sets  the  count  to 
zero,  if  two  successive  terras  (s  is  the  value  of  the  preceding  term) 
are  of  the  same  sign.  When  the  count  succeeds  in  reaching  a  value 
greater  than  the  Integer  i(l),  also  used  in  Step  1.44*,  the  process 
is  presumed  to  have  shown  that  further  summing  of  the  series  is 
unnecessary  and  control  passes  to  Step  1  .<>»>.  If,  however,  the  count 
does  not  exceed  1(1),  Step  1.54  send  »  Che  control  to  Steps  1.58 
through  1 . o4  in  which  the  first  three  stops  are  partial  preparation 
for  the  next  roet  in'*  le  in  the  r  >w  bet  ore  Step  1  .84  passes  the  control 
t>  Step  1.12,  where  the  proper  »t  I  on  is  .-■>  -n  le  tod  . 


iihvi  mix  !  v ,  it  experience  in, lie. lies  ditteriug  values  should  he  used,  the 
i  ;  no  r* ■  isnti  whv  n  >(  . 
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The  Algorithm  tor _ the  cos  2kt_-  eos_  2_kq-  Component  of  the  Second  Integral 

Because  cos  2kx-  vanishes  for  x  =  fi  /  2 ,  3't/2,  5n/2,  . ..,  the  width  of 
the  first  row  and  also  the  first  column  is  [/(4k)  and  the  lines  subdividing 
the  area  of  integration  are 

t  =  V (2n-l )*p 

and 

q  =  V (2n-l)*p  , 

where  p  =  t/(4*k),  (2n-l)  =  1,3,5,...  and  (2m-l)  =  1,3,5,...  . 

>  2 

Only  a  few  steps  in  the  algorithm  tor  the  sin  2kt  "  sin  skq“  component 


chair 

The  steps 

are 

*, 

1.1 

Set 

P 

■t  /  C  4 

*k) 

1.26 

Set 

3(2)  = 

sort 

l(2*N(2)  - 

l)*pl 

1.32 

Set 

BCD  = 

sqrt 

l ( 2*N(  1 )  - 

1) *p] 

The  ,\L  torithm  tot  the  sin  2kt~  cos _ 2kq-  Component  of  the  Se_cond_  LnteiuraJ. 

For  this  component  the  width  of  the  first  column  is  r/(2k)  and  of  the 
first  row  t/(4k.).  Two  p-values  are  needed 

PC 1 )  =  */(2*k) 

and 

p(2)  =  t/(4*k)  . 


Hence 

l  .  1 

Set 

p(l) 

=  t/(2*  k) 

1.12 

Set 

p(  2 ) 

=  -t/  (  4*k) 

1  .26 

Set 

b(2) 

=  sqrt  l(2*N(2)  -  l)*p( 2) ) 

1.32 

Set 

BCD 

=  sqrt  [N(  1  )*p( 1 ) ] 

I 
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The  Algorithm  for  the  cos  2kt^  sin  2_kq-  Component^  ot  the  Second  Integral 
In  this  case 

p(  t  )  =  -i/(4k)  , 


and 


p(2)  =  :i/(2k)  , 


so 

1 . 1 

Set 

PC  l ) 

1.12 

Set 

p(  2 ) 

1  .  2b 

Set 

B(2) 

1.32 

Set 

B(l) 

The  Algorithms  for  the  First _ Integral 

The  first  integral  differs  from  the  second  in  its  area  of  integration 
and  its  integrand.  The  area  ot  integration  of  the  former  is  a  strip  in  the 
t-direction  with  a  finite  boundary  at  t  =  l.  The  area  of  integration  of  the 
latter  is  the  same  except  that  the  finite  boundary  is  at  t  =  o.  This 
difference  is  of  no  great  consequence. 

The  difference  in  the  integrands  is  a  little  more  significant.  The 
first  integral  has  the  quant i ties  V  t  - 1  and  q  -  t  in  its  integrand,  the 
second,  the  quant! Cites  \t+ l  and  q  +  t. 

The  infinite  singularity  originating  in  v  t  -  l"  is  integrable  and 
readily  handled  by  means  ot  quadrature  formulas  for  the  weight  function 

r 

i  r\t. 

The  singularity  caused  by  q  -  t  is  not  integrable,  but  it  appears  in 
tile  corner  of  the  area  of  integration,  where  q  =  t  =  1  and  where  the 
quant i t v  y  q( 1  -  q )  vanishes.  The  improper  double  integral  can  be  shown 
to  exist  by  comparing  it  with  the  improper  double  integral  obtained  bv 
replacing  the  trigonometric  factor  by  unity.  The  singularity  is  no  problem 


numerically,  though  it  does  slow  the  evaluation. 
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Tile  changes 


to  convert  algorithms  tor  the  second  integral 


into  those  for  the  first  are  few.  In  the  case  where  the  trigonometric 
n  ■) 

factor  is  sin  2k.t-  sin  2kq“  thev  are: 

1.12  Set  n  =  first  [i  =  l(l)l't)  ;  i*>>  >  1] 

1.2  Set  N( L )  =  n 
1.22  Set  A(l)  =  l 

2.31  To  Step  2.5'3  if  1  =  1  and  i'lI,U(.l)|  =  1 
2.325  To  Step  2.34 

2.51  Set  X(l)  =  I'll, 0(1)]  +  d  [  J  (  I  )  ,  K(  I  )  ]  *  H 1  I  ,  1  ( 1 )  ] 

2.33  To  Step  2.57  if  l  =  1  and  H[  l, <>(!)]  =  1 
2. 5b 5  To  Stop  2.58 

2.57  Set  V[  i  ,.J(  I)  ]  =  V  [  I ,  J  (  1  )  ]  =  c  |  .1  (  1  ) ,  l  )  ]  *  sqrt  >[1,0(1)]  *  G(I) 

The  right-hand  member  of  Step  1.12  is  a  function  or  i  that  finds  the 
first  of  a  .sequence  of  values  of  i  that  neets  the  prescribed  condition. 

In  this  case,  the  sequence  is  i  =  1 , 2  ,  3 ,  . .  .  , 100  and  the  prescribed  condition 

is  that  p  *  i  be  creator  than  unity.  The  value  of  n  for  each  value  of  k 

can,  of  course,  be  found  by  a  brief  calculation  and  treated  as  an  input. 

One  can,  however,  elect  to  write  a  small  subroutine  equivalent  to  the 
function  of  Step  1.12. 

The  steps  added  to  Part  2,  it  wi L 1  be  seen,  provide  an  alternative  to 
Gaussian  quadratures  whenever  I  =  1  and  T( 1,0(1))  =  1.  The  alternative  is 


a  Gaussian  quadrature  for  the  weight  function  1 


!"\jt.  The  lower  cast 


symbols,  c  and  J,  are  respectively  the  weights  and  abscissas  corresponding 
to  the  latter  weight  function  on  the  interval  [0,1].  They  are  derivable 
from  the  weight  and  abscissas  for  w(x)  =  1  in  the  interval  (-1,1).  See 
Appendix  B. 

The  changes  necessary  to  convert  the  remaining  three  algorithms  for 


tlie  second  integral  into  those  for  the  first  integral  differ  very  littb 


from  those  above. 
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The  steps  added  to  Part  2  are  stiLl  needed  to  take  care  or  the 

singularity  at  t  =  1  and  must,  of  course,  be  retained.  The  only  step  that 

2  > 

changes  in  Part  1  is  Step  1.12.  For  the  factors  sin  2kt~,  cos  2kq“,  one 

has  1.12  n  =  first  [i  =  1(1)100:  i*p(l)  >  1].  The  only  change  is  the 

*>  •) 

unit  index,  or  subscript,  for  p.  For  the  factors  cos  2kt~,  cos  2kq"  and 

i  > 

cos  2kt~,  sin  2kq~  the  step  defining  n  becomes 
1.12  n  =  first  li  =  1(1)100:  (2*i  -  l)*p(l)  >  1]. 

Step  1.12  now  reads  n  =  the  first  value  of  i  in  the  sequence 

i  =  1 , 2 ,),...,  100  such  that  (2*i  -  l)*p(l)  >  1  anti  n  =  1,5,3 . 199. 

The  complete  FORTRAN  program  for  computing  the  lag  function  is 
presented  in  Appendix  C. 
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Appendix  A;  FORTRAN  Algorithms 

The  following  algorithms  are  the  FORTRAN  coded  equivalents  of  the  first 
three  algorithms  appearing  in  the  text.  Statement  numbers  occurring  in  the 
text  have  been  retained,  where  possible,  to  facilitate  easy  following  of  the 
logical  path  of  the  programs. 

Algorithm  3  has  been  written  as  two  subroutines.  The  first  performs 
the  "outer"  iteration  which  depends  on  the  value  returned  by  the  second  or 
"inner"  iteration.  The  subroutines  would,  in  practice,  be  called  by  a  main 
program  which  passes  the  values  of  the  "outer"  limits  of  integration  [A,B] , 
the  "inner"  limits  [A2,82]  and  a  parameter  R  (if  needed  for  the  evaluation  of 
F(X,Y,R)j.  The  result  of  the  double  integration  is  the  variable  S  returned 
to  the  main  program. 


-29- 


27  March  1981 
GHP : SJG : Ihz 


1 


ALGORITHM  1:  2,  3,  4  POINT  GAUSSIAN  QUADRATURE  ON  FULL  INTERVAL  [A,  B] 


PROGRAM  GAUSS  I 
IMPLICIT  REAL-3 ( A-HtO-Z ) 

DIMENSION  C  (4»4)  »D  (4,4)  ,V  (4-) 

DATA  E/.0002/, A/0.0 / t 3/3. 141 5926 5/ TC( 2. l)/.500/,D( 3,2 )/.500/ 
DATA  C(3,l)/.277777777777777D0/TD(4»l)/.06943l84420297300/ 
DATA  C{  3»2) /.444444444444444DO/  ,C  (4,2  )/. 3 2 50 72 5 774 3 1 2 7 3DO/ 
DATA  C(4» 1) /. 17392742256372700/ ,0 ( 3tl)/.ll2701665379253D0/ 
DATA  0 ( 4t 2 )/. 33 000947S 20757200/, I  PR/3/ 

F(X)=DSIN(X) 

C{2t2)=C(2t1) 

C(3t3)=C(3tI) 

C (4,3)=C(4t 2) 

C(4»4)=C{4,1) 

D(2tL)  =  . 5D0  — CSG  RT ( 3  >000 ) /6.0D0 
D(2,2)=1-ODO-Q'(2tI) 

D{3,3)=1.000-D{ 3tL) 

0 (4»3 )=1.00Q-D(4t2 ) 

D(4t4)=1.000-0(4t1) 

H  =  9— A 

V(1)=0-0D0 
DO  5  J  =  2 » 4 
V{ J)=0-000 
DO  7  K=1tJ 
X  =  A«-D(  J,K)AH 

7  V{ J)=V{ J) ♦C ( JtK)-HAF(X) 

IF (DABS (V ( J )— V ( J-l ) ) .  L  E  -  E  )  GO  TO  13 

5  CONTINUE 

13  WRITE(IPRtI)  J»V(J) 

1  FORMAT! I2»2X,E16-3) 

STOP 

END 


i 
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ALGORITHM  2:  2,  3,  4  POINT  GAUSSIAN  QUADRATURE  WITH  SUBDIVISION  OF  INTERVAL  [A,  K] 


PROGRAM  GAUSS2 
IMPLICIT  REAL* 9 ( A-H,0-Z ) 

DIMENSION  C(4,4)»D(4»4)»V(4)»H(10)»U(10) 

DATA  E/. 0C02/, A/O. 0/, 3/3. 14 1 59265/, C( 2,1)/. 5D0/,D( 3, 2)/. 500/ 
DATA  C(3,1)/.2T77T7777777777DO/,0(4,  1) /.0694 3  I  344202 97 3 00/ 
DATA  C { 3, 2)/. 44444444444444400/, C {4,2 )/. 3 260 72 5 7 74 3 L 27 3 00/ 
DATA  C (4, 1)/. 17392 742 2 56 8 7270 0/,0{ 3 ,1)/.112701665379253D0/ 
DATA  D(4,2}/-3300094782Q7572D0/,IPR/3/ 

F ( X )  =DS  IN ( X ) 

C(2,2)=C(2,1) 

C (3,3)=C( 3, 1) 

C (4,3)=C(4,2) 

C ( 4 , 4 ) =C ( 4 , 1 ) 

0(2,1)  =  .500-0 SORT  (3. 000)'/6.000 
D{2,2)=1. 000-0(2,1) 

D(3,3)=l. ODO-D (3,1) 

0(4, 3)=1. 000-0(4, 2) 

D(4,4)=l. 000-0(4,1) 

N=1 

H ( N )  =  B— A 
U  (  N  )  =  A 
5=0.000 

5  V(1)=0.000 

DO  7  J  =  2 , 4 
V( J)=0-000 
DO  10  K  =  1 , J 
X  =  U  ( N  )  ♦  0  (  J  ,  K  )  *H  ( N ) 

10  V ( J )=V{ J) *C ( J ,K )*H(N)*F (X ) 

IF (DABS (V (J)-V(J-l ) }«LE.E)  GO  TO  21 
7  CONTINUE 

H(N)=H(N )/2.0 
N  =  N*1 

H(N)=H(N-1 ) 

U(N)=U(N-1)*H(N) 

GO  TO  5 
21  S=S*V(J) 

IF (N.EO.l )  GO  TO  25 

N  =  N—  1 
GO  TO  5 

25  WR  I  TE { I  PR  ,  l  )  J,S 

1  FORMAT( 12, 2X, 616. 8) 

STOP 

END 
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ALGORITHM  3:  TWO  DIMENSIONAL  2.  3,  A  POINT  GAUSSIAN  QUADRATURE  WITH 
SUBDIVISION  OF  INTERVAL  [ A , B ] 


SUBROUTINE  GAUSSA( A,B,A2,62,R,S) 

IMPLICIT  REAL-8 ( A-H,0-Z ) 

DIMENSION  C  (  4  » 4  ) »  D  (  4 , 4  ) ,  V (  4  )  ,  H  { 10) ,  U  ( 10) 

DATA  E/. 0  00 200/ ,C (  2 ,  1 ) / . 5 DO/ , D ( 3 , 2 ) / . 500/ 

DATA  C(3,1)/.277777777777777D0/, 0(4, l)/„0694 3184420297300/ 
DATA  C( 3, 2 )/.444444444444444D0/,C ( 4 , 2 ) / . 3260 72 5774 3 1 2 7 3D0/ 
DATA  C  { 4  ,  1 )/.  17  392  74  22  568  72700/,  0(  3, 1 )/.  1 1  270  166  53  792  58  00/ 
DATA  0(4, 2)/. 33 000 947 8 207 57200/, I  PR/ 3/ 

F ( X,Y,R )=Y-DS IN ( X ) 

C(2,2)=C(2,1) 

C(3»3)=C(3,1) 

C (4,3)=C(4,2) 

C (4,4)=C(4,1) 

0(2,1  }=. 500-0 SORT ( 3 .000 ) / 6. ODO 
D(2,2)=1.0D0-D(2,1) 

0(3, 3)=1. 000-0(3,1  ) 

D(4,3)=l. 000-0(4, 2) 

0(4,4)=1.000-D(4,l) 

N  =  1 

H ( N)=B-A 
U  (  N  )  =  A 
S=0-0D0 

118  V ( 1 ) =0.000 

DO  138  J  =  2 ,4 
V( J)=0.0D0 
DO  150  K=1»J 
X=U(N)+D(J,X) -H ( N ) 

CALL  G1(X,Y,A2,B2,R) 

150  V(J)=V(J)+C(J,K)-H(N)-F(X,Y,R) 

IF (DABS (V ( J)-V( J-l ) ) .LE.E)  GO  TO  134 
138  CONTINUE 

H(N)=H(N)/2.0DO 
N  =  N+  1 

H  (  N  )  =  H  (  N—  1 ) 

U  (  N ) =  U ( N— 1 ) +  H  (  N  ) 

GO  TO  118 
184  S  =  S ♦V ( J  ) 

IF(N.EO.l)  GO  TO  194 
N  =  M-1 
GO  TO  118 
194  RETURN 
ENO 


-32- 


27  March  19bl 
GHP:SJC:lhz 


ALGORITHM  3:  TWO  DIMENSIONAL  2,3,4  POINT  GAUSSIAN  QUADRATURE  KITH 
SUBDIVISION  OF  INTERVAL  [A,B]  -  continuation 


SUBROUTINE  G I ( Y » S , A , B,R ) 

IMPLICIT  REALMS (A-H,0-Z ) 

DIMENSION  C(4,4),D(4,4),V(4),H(10),U(10) 

DATA  E/ . OOO  200/ ,C ( 2 , 1 ) / . 500/ , D ( 3 t 2 ) / . 500/ 

DATA  C( 3* 1  )/.27777777777777700/,D{4, 1 ) /. 06943 1S44 202 973D0/ 
DATA  C ( 3,2 )/. 44444444444444400/ ,C (4,2 )/. 3260725774 31 27300/ 
DATA  C (4,  1 )/- 17 3927422 568727D0/,D ( 3, 1  ) /. 1 1 27 0 1 66 5 3 792 5800/ 
DATA  D (4, 2)/. 33000947820757200/, I  PR /3/ 

F(X,Y,R)=DSIN(X) 

C{2,2)=C(2,1) 

C ( 3  » 3 ) =C ( 3 , 1 ) 

C(4,3)=C(4,2) 

C(4,4)=C(4,1) 

D  (2, 1 )=.500-D SORT (3. 000 )/6.000 
D(2,2)=l. 000-0(2,1) 

D(3,3)  =  1.0D0-D( 3,  1  ) 

D (4,3 }= 1.000-0(4, 2 ) 

D (4,4 )  =  1«0D0-D(4, 1  ) 

N=1 

H(N)=B-A 
U  (  N  )  =  A 
5=0-000 

118  V(l)=0-000 

DO  138  J=2,4 
V ( J )=0-0D0 
DO  150  K=  1 , J 
X=U(N)*0( J,K)*H(N) 

150  V( J)=V( J) *C ( J,K)*H(N)-F(X,Y,R) 

IF (OASS(V(J)-V(J— 1 ) ).LE.E)  GO  TO  184 
138  CONTINUE 

H(N)=H(N)/2-0D0 
N  =  N*  1 

H(N)=H(N-1) 

U  (  N  )  =  U  (  N—  1  )  ♦  H  ( N  ) 

GO  TO  118 
134  S=S*V(J) 

IF(N.EO.l)  GO  TO  194 
N  =  N-  1 
GO  TO  118 
194  RETURN 
END 


I 
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Appendix  13 :  Two  Sots  of  Gaussian  Weights  and  Absci  ssas 


INTERVAL  =  [0,1] 


■h 


.5 

.5 

n=3 

.27777  77777  77778 
.44444  44444  44444 
C3  =  CL 

n=4 

.17392  74225  68727 
.32607  25774  31273 
C3  =  C2 
C4  =  Ci 

n=5 

.11846  34425  28095 
.23931  43352  49683 
.28444  44444  44444 
C4  =  C’2 
C5  =  Ci 

n=6 

.08566  22461  89585 
.18038  07865  24070 
.23395  69672  86345 
C4  =  c3 
C5  =  C2 
C6  =  Cl 


WEIGHT  FUNCTION  =  1 


;<1 


.21 132  48654  05187 
.78867  51345  94813 


.11270  16653  79258 

.  5 

.88729  83346  20742 


.06943  84420  29730 
.33000  94782  07572 
.66999  05217  92428 
.93056  81557  97027 


.04691  00770  30668 
.23076  53449  47159 

.5 

.76923  46550  52341 
.95308  99229  69332 


.03376  52428  98424 
.16939  53067  66867 
.38069  04069  58401 
.61930  95930  41599 
.83060  46932  33133 
.96623  47571  01576 
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INTERVAL  =  [0,1]  WEIGHT  FUNCTLON  =  1 /yjH 


n=2 

11558 

71099 

97048 

1.3042 

90309 

i  72509 

74155 

57471 

45810 

.69570  96902 

:  74908 

n=3 

05693 

11596 

70074 

.93582 

78691 

45382 

43719 

78527 

51095 

.72152 

31460 

96278 

86908 

43784 

32472 

.34264 

89847 

58340 

n=4 

03364 

82680 

67507 

.72536 

75667 

56724 

27618 

43138 

72464 

.62741 

32917 

55774 

63467 

74762 

34637 

.44476 

20689 

06748 

92215 

66084 

92058 

.20245 

70725 

80752 

n=5 

02216 

35688 

07218 

.59104 

84494 

29951 

18783 

1  5676 

52445 

.53853 

34386 

19992 

46159 

73614 

96266 

.43817 

27250 

31964 

74833 

46283 

87281 

.29890 

26983 

01162 

94848 

39262 

28836 

.13334 

26886 

17376 

n=6 

0 1  568 

34066 

07400 

.49829 

40916 

26806 

13530 

00116 

55248 

.46698 

507  30 

76710 

38069 

04069 

58401 

.40633 

48534 

46132 

61930 

95930 

41599 

.32013 

66570 

86692 

81742 

80132 

66875 

.21387 

36519 

90636 

96346 

1296  3 

70913 

.09435 

06727 

7  3024 
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Appendix  C:  implementation  of  Algorithms  for  Computing  the  Lag  Function 

The  foLLowing  example  is  an  adaptation  of  the  preceding  algorithms 
applied  to  the  computation  of  the  lag  function.  Implementation  follows 
directly  from  that  outlined  in  the  main  body  of  this  manuscript.  The  few 
exceptions  are  clarified  in  the  following  outline. 

1.  MAIM  PROGRAM :  For  each  value  of  R  (R=k),  eight  double  integrals  are 
evaluated.  The  real  and  imaginary  parts  (SRF,  SIM  respectively)  are  the  suns 
of  the  corresponding  contributions  (SUM(l)  )  of  the  appropriate  Integrals.  No 
output  statements  are  supplied  as  the  relevant  data  are  subject  to  user 
discretion. 

For  efficiency  and  economy,  the  rj  interval  (0,1 j  was  held  constant  in  the 
main  program  since  the  q  algorithm  provides  for  subdivision  of  the  interval 
when  necessary.  In  light  of  this,  a  six-point  formula  was  used  which  is 
theoretically  exact  for  polynomials  of  order  eleven  or  less.  For  the  maximum 
value  of  k  used  (10.0),  the  number  of  nodes  on  the  interval  [0,1]  is  less  than 
eleven. 

T(l)  and  T(2)  are  the  n  and  (n-1)  partial  suns  of  the  alternating  series . 
'Then  the  series  is  terminated,  they  provide  an  upper  and  lower  bound  on  the 
solution.  T(3)  is  the  midpoint  and  is  accepted  as  the  correct  result.  T(4) 
and  T(5)  are  summations  of  the  total  estimated  errors  in  the  t  and  u 
integration  respectively. 

E(1MT)  is  used  as  a  pseudo  relative  error  bound  on  the  n1'1  partial  sum. 
When  the  contribution  to  T(l)  is  less  than  K( LNT)  or  some  predetermined  value 
(  .0001  in  this  case),  the  summation  of  the  alternating  series  terminates. 
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l l .  SUBROUTINES: 

(a)  GAUSS1 ,  GAUSS2  -  These  perform  the  quadratures  for  t  and  q  respectively. 
The  weights  and  abscissas  for  each  type  of  quadrature  are  initialized  in 
GAUSS  1  and  passed  to  GAUSS2  in  the  argument  List  of  the  CALI,  statement.  Thes 
subroutines  are  the  FORTRAN  equivalents  of  the  last  algorithm  in  the  text 
which  uses  a  relative  error  bound.  Variables  not  appearing  in  the  algorithm 
are  explained  as  follows. 

.JJ(N)  saves  the  value  of  the  order  of  the  last  quadrature  accepted.  Thi 
prevents  duplication  of  effort  in  recalculating  lower  order  quadratures 
previously  computed  and  accepted  by  GAUSS  1  and  GAUSS2. 

SV(N)  saves  the  contribution  of  the  N*-'1  subinterval  which  is  added  to  th 
total  sum  at  Step  100  should  the  condition  at  Step  200  be  met. 

F.Q1  is  the  estimated  error  in  the  q-integration  for  a  specific  value  of 
t.  TEQ  is  the  sum  of  the  estimated  q  error  in  the  entire  t  interval  just 
computed . 

(b)  PVAL  -  Calculates  the  factor  used  to  determine  the  zeroes  of  the  sin  or 
cos  argument  for  the  boundaries  of  the  subintervals  on  the  infinite  t  range. 

(c)  LVAL  -  Determines  the  integer  multiplier  for  the  first  upper  limit  in  t 
when  the  lower  limit  =  1.0. 

(d)  AIV  -  Sets  the  lower  limit  on  the  t  interval  to  0.0  or  1.0  depending  on 
which  integral  is  being  evaluated. 

(e)  BIV  -  Calculates  the  next  node  point  on  the  t  interval. 

(f)  FT/FQ  -  Evaluates  the  function  at  a  specific  abscissa  in  the  appropri  ite 
t/q  quadrature. 


120 


132 


156 

156 

138 

900 

800 


main  program  ***0* 

IMPLICIT  REAL1'?  (A-H.Q-Z  ) 

0  I  MANSION  T  (  5  )  .  SON  (  3  )  ,E  (  6  )  =  (  34  ) 

COMMON  AI  ,3  I  ,  A2  ,52  .E  *  If-'T 

DATA  1 1/5/ . RE/. 0005  2)9/, I »?/&/.? I / i.  IM  5M 
DATA  RR/.Gl,.02,.04,.06,.03,.l,.12,.la.. 

>8. 0,9. 0,10.0/ 

DO  800  NR=l,24 
P -RR ( NR ) 

DO  900  I N  T  =  1 , 8 
SUM ( I  NT )  =  0  *  ODO 
S=0.000 
A  2  =  0  • 000 

CALL  PVAL ( PI , P2 »P  I  ) 

CALL  L V  AL ( P 1 »  L ) 

I  1  =  L  +  4 
M  =  l 

CALL  A  1  V 
T (  1  ) =0.0D0 
T (4) =0.000 
T  (5)  =  0.000 
82=1.000 
CALL  B1V{M,P1) 

CALL  GAUSS1  ( SI, 6T,EQ) 

T ( 4 ) =T ( 4 ) ♦ET 
T (5)=T( 5) -EQ 
T(2)=T( 1) 

T(1)=T(1)*S1 

E ( INT)=OMAX 1 ( 0A3S (T( 1 )*RE ) » .10-3 ) 
IF(M.LT.Il)  GO  TO  156 
IF (DABS (SI ) .GT. E{ INT) )  GO  TO  156 
IF(S*S1.GT. 0.000)  GO  TO  156 
IZ=IZ*1 

IF(IZ-Il)  158,158,188 
I  Z  =0 
S  =  S1 
A  1=31 
M  =  M  *  I 
GO  TO  132 

T(3)=(T(1)*T(2)  ) / 2 .000 
E ( INT )=0A3S ( SI) /2.000 
SUM ( I NT )  =  T ( 3 ) 

CONTINUE 

SRE  =  SUM ( 1 ) *  SUM ( 2 ) *  SUM (3) ♦SUM (4) 

S I  Ms SUM ( 5 ) -  SUM ( 6 ) *  SUM ( 7 ) -  SUM ( 3 ) 

CONTINUE 

STOP 

ENO 
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SUB ROUTINE  G  a  u  S  S  1 (S»TE»TEG) 

I  °  L  I  C  I  T  REAL-3 ( A-H,0-Z ) 

COMMON  A, B , A2 ,32 , R  ,  I NT 

DIMENSION  C  {  6 ,6 ) , D ( o  ,6  )  ,7 ( 6  )  ,  H  (  1  50  )  ,  U  (  l 50 ) , 0 (  1  50  )  ,  W ( 6 , 6  )  , Z ( 6 . 6  ) 

01 "  =  NS  I  ON  SS(  150)  ,  J  J  (  15C)  ,  S  V  (  150) 

DATA  C  (  2  ■»  1)  ,  C  (  2  *  2  )  /  2  .  5  0  0  / 

DATA  C  (  3  *  1)  *  C  (  3  »  3  ) ,  C  (  3  »  2  ) 

>/ 2  ■■■.Z  77777777777777300,  . 44 A44A444444 444400/ 

DATA  C  {  4 , 1  )  »  C  (  4  »  4  )  »  C  (  4 , 2  )  »  C  (  4  ,  ?  ) 

>/ 2-. 17392 74 2 2 56 372 7 00,  2 *. 3260 7 2 57743 1 2 7 300/ 

DATA  C(5,1),C(5»5) *  C ( 5  »  2 ) *C(5*4)»C{5»3) 

>/2*.l  18  46344  2  52  30  95  2  CO,  2- .  2  3^  3  1  4  3  3  5  2  49  6  3  3D0  , 

>  .23444444444444400/ 

DATA  C(6,L),C(6,6),C(6,2)»C(6»5),C(6»3),C(6,4) 

>/2-.  356622461395350-1,  2:-:.  1303807  3652407000,  2  *.  2  3  39  56s  6  7  2  36  3  4  5  00/ 
DATA  0(2,1)  ,0(2,2) 

>/ . 2 l 1324865405137100 ,  .789673134594312900/ 

DATA  0(3,1), 0(3, 2),  0(3, 3) 

>/. 11270166537925900,  .500,  . 3 6 7 29 3 3 3 46 20 74 2 DO/ 

DATA  D(4,  1)  ,0(4,2)  ,0(4,3)  ,0(4,4) 

>/.o94 3 134420297 30- 1 ,  .330009473  20  7  572  00,  .66  99  9052  179  242300, 

>  .93056315579702700/ 

DATA  0(5,1), 0(5, 2), 0(5, 3), 0(5, 4), 0(5, 5) 

>/ .469100770306630- 1  ,  .2307653449^715900,  .500, 

>  .76923465505284100,  .  9 5 30 89 9 2 29o9 3 3 2 00/ 

OATA  0(6,1), 0(6, 2), 0(6, 3), 0(6, A), 0(6, 5), 0(6, 6) 

>/. 33765 24 29934240-1  ,  .  169  39  5  30  67  66  367  00,  .  33069Q 40 695 3 40  1  DO , 

>  .61930959304159900,  .83060449323313300,  .96623475710157600/ 

DATA  W ( 2.  1  )  ,W ( 2 ,2 ) 

>/ 1.304290 3097 2509 2 00  ,  .695  70  96  90  27490300/ 

DATA  W( 3, 1)  ,W(3,2) ,W( 3,3) 

>/. 93532736914533200,  . 7 2  1  5 2 3  1 4609 6 2 7 8 00 ,  .342643  93^7  5  8  3  400  0/ 

DATA  W ( 4 , 1 )  ,  W ( 4 , 2 )  »  W ( 4  »  3 )  ,  W  ( 4 , 4 ) 

>/. 72536756675672400,  .62741329175577400,  .44476206590674300, 

>  .20245707253075200/ 

DATA  W ( 5 , 1  ) , W { 5 , 2 ) , W ( 5 , 3 ) , W ( 5 , 4 ) , W ( 5 , 5 ) 

>/. 591 04 34 49 42 99 50600,  .53353343861999200,  .433 17272503 196400, 

>  .29390269330116200,  .13334263361737600/ 

OATA  W ( 6 , 1  ) ,W(6,2)  ,4(6,3)  »  W  (  6  *  4  )  ,4(6,5)  ,  W  (  6 , 6 ) 

>/. 49829409162630600-  .46693507307671000,  .40633435344613200, 

>  .320  15665703669200,  .2  1  337865199063600,  .  943 50 6 7 2 77 3 0 2 40-  1  / 

DATA  Z ( 2 , 1 )  ,  Z (2 ,2 ) 

>/.l 15537109997047700,  .741555747145309900/ 

OATA  Z ( 3 , 1 )  ,  Z ( 3 , 2 ) , Z ( 3 , 3 ) 

>/. 56939  11 59 67 00 740-1,  .437  1973527510946  00,  .86^03437343247  1500/ 

DATA  Z(4»l) »Z(4,2) , Z ( 4 , 3 ) , Z (4 ,4 ) 

>/. 3 36 43 2 6306 7 50 69 2 0-1,  .276134313872464400,  .634677476234637100, 

>  .922  156603492053300/ 

OATA  2(5,1) ,Z (5,2) ,Z(5,3) ,Z(5»4) ,Z(5,5) 

>/. 22163 56  33072  1757  0-1,  .  1370  31  56  7652445300, 

>  .461597361496266200,  .748334623337231300,  .943493926238369100/ 

OATA  Z(6,1),Z(6,2),Z(6,3),Z(6,4),Z(6,5),Z(6,6) 

>/. 156334066074004 70-1,  .135300011655248100,  .33069040695340100, 

>  .61930959304159900,  .317423013266375300,  .963461294370912800/ 

Data  RE/. 50-3/, V( 1 )/0. 0000/ 
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CHP :  SJC. 

s  =  1 
M=  1 

jj  ( n  =2 

SV (  1  ) =0.0 CO 
0  {  N)  =  I  .00t>3 
E  =  l  .0062 
H  (  M)  =3-4 
IJ  {  M)  =  A 

2  22  S  =  0 .000 

IE =0. GOO 
TE'Z-O.ODO 

2  DO  IF{0(N)-lE.E)  GO  TO  100 

N  J  =  J  J  (  N  ) 

Ir(.'JJ-o)  30,30,500 
30  DO  10  J  =  N  J  »  6 
V ( J ) =0.000 
EQ1=0.000 
DO  20  K-ltJ 
X=U(N)+0( J,K)-H(N) 

GO  10(1,1,2,2,1,1,2,2)  ,[K'T 

1  I  -  {  U(  N)  .EQ.  1 .000)  X=U(N)  +Z(  J,K)*H(M) 

2  CALL  GAUSS2(X,Y,C,0,EQ) 

E  0 1  =  E  0 1 +-EQ 

GO  TO (3, 3, 4, 4, 3, 3, 4, 4) ,  INT 

3  IP (U(N) .EO. L.OOO)  GO  TO  400 

4  CALL  FT (  X  ,  Y  ,  F  X  ) 

V ( J ) =  V ( J ) *C ( J , K ) *H ( N ) *FX 
GO  TO  20 

*-00  CALL  FT  (X,Y,GX) 

V  (  J  )  =  V  (  J  )  ♦  W  (  J  ,  K  )  0  5  0  ?  T  (  2  (  J  ,  K  )  ) :':  H  (  0  )  *  Q  X 
20  CONTINUE 

S  S  (  N  )  =  V  (  J  ) 

J  J  (  N  )  =  J  *  1 

Q ( N) =0A35 ( SS ( N ) - S V ( \ )  ) 

IP(0(N) .LE.E)  GO  TO  100 
S  V  (  N )  =  S  S  {  N ) 

10  CONTINUE 

500  H(N)=H(N)/2.000 

W  =  M  *  1 
H(M)=H(N) 

U( M)=U(N) 

Q(N)=Q{N) 

JJ (N) =2 
JJ(*)=2 
SV(N)=0.0D0 
S  V ( M ) =0.000 
GO  TO  200 
100  S  =  S  *  S  S  (  N ) 

SV ( N) =SS ( N) 

Tc=TE*0(N) 

TE0=TE0*E01 
N  =  N*  1 

IF (N.LT.H-l )  go  TO  200 
IF(TE.LE.RE-OABS(S) )  GO  TO  300 
N  =  l 

£=  (TE-RE-0A2S  (  S  )  ) /,M 
GO  TO  222 
300  RETURN 

END 
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Sl'^.O'OT  [M=  GAUS32(''.S.C,0.T6) 

IMPLICIT  (  A  -  u  ♦  0  -  7  ) 

C  3  -‘NON  A  i  ,  9  1  ♦  A  .  3  .  i  .  2  T 

Di,J:':S  ION  C  (6  ,6  )  ,  3  (o  ,5  )  ,  V  (6  )  ,H(  1  50  )  ,'J  (  150  }  ,0  (  150) 
DIMENSION  S3 (  150)  t JJ (  150)  tSV (  150) 

0  A  7  A  RE/.5D-4/* V( i )/0.DD00/ 

N=l 

JJ ( 1 ) =2 
SV (  1 ) =  0 . 0  DO 
£= l .0052 
0( N) = 1 .0063 
”  =  1 

F('l)=8-A 
U  (  N )  =  A 

222  S=0.000 

TE=0.000 

200  IF(Q(N) .LE.E)  GO  TO  100 
NJ  =  JJ (N) 

IF(NJ-6)  30»30,500 
30  DO  10  J  =  N  J  *  6 
V { J ) =0.000 
00  20  K=1,J 
X  =  U  (  N  )  *■  D  (  J  ■*  K  )  >:  H  (  N  ) 

CALL  FQ(X,Y.FX) 

2  0  V(  J)  =  V(  J)  *C  (  J*FX 

SS(N)=V{ J) 

J  J  (  M  )  =  J  *  l 

Q  ( N )  =  0  A  3  S  (SS(N)-SV(N)  ) 

IF (0(N)  .LE.E)  GO  TO  100 
SV(NJ sSS(N) 

10  CONTINUE 
500  H(N)=H( N)/2.000 

M=M*  1 
H(M)=H(N) 

U{M)=U(N) ♦  H  (  N  ) 

0(M)sQ(N) 

J  J  (  N  )  =  2 
J J ( M) =2 
SV(N) =0.000 
SV(M) =0.000 
GO  TO  200 
100  S  =  5  *  5  S ( N ) 

S  V  (  N ) =SS(N) 

TE=TE*Q{ N) 

N  =  N*  l 

IF  (N.LT.M*!)  GO  TO  200 

Ir{ Tc.L5.RS*0A8S(S) )  GO  TO  300 

Ns  1 

E= (T5-RE-QA3S (S) )/M 
GO  TO  222 
200  RETURN 
END 
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1 

2 

3 

4 


SU3R0UT  IN 

E  p  y  a  l 

(  =  ] 

L  ,  P2  , 

°l ) 

I MPLICI  T 

REAL-0 

(A- 

•M  ,  0- 

? ) 

CO^'ON  A  1 

,  3  1  .  A  2 

*  ^  i 

?  »  R  ,  I 

\T 

GO  TO (1,2 

,1,2,3 

•  3,4) 

,  I  N  T 

°l=PI/( 2. 

OOC-R) 

P2=PI/(2. 

CDD-R) 

R ETUR N 

p 1=P I / ( 4.CD0-R ) 
P2=°I/{4.0DC4R) 
RETURN 

P1=PI/(4.CD0*R) 
P2=PI/{ 2.0DC-R) 
RETURN 

Pl=°I/(2-  ODO*R  ) 
P2  =  3I/(4.00C-::R) 
RETURN 
END 


SUBROUTINE  LVAL { P 1 ,L  ) 

IMPLICIT  REALMS (A-H,0-Z) 

COMMON  A1,31,A2,52,R,INT 
GO  TO ( 1 ,2,3,3 ,2 ♦ 1 ,3,3 ) - 1  NT 

1  DO  10  1=1,100 
L  =  I 

if(l*pi.gt.i.cqo)RETurn 

10  CONTINUE 

RETURN 

2  -CO  20  1=1,100 

L  =  I 

I  F  ( (  2  -  L  - 1 )*P1«GT. 1.000)  RETURN 
20  CONTINUE 

RETURN 

3  L  =  1 
RETURN 
END 


<;■>  n 


C! 


SUBROUTINE  A  1  V 
ILLICIT  REAL*;  (A-H, 0-2) 
0M.*ON  Al  »31  .  A2  ,E2  ,R  ,  I  mt 
G  TO  (2 ,2  ,  1  ,  1  ,  2 ,2  ,  1  ,  1  )  ,  I  NT 
ai=o.oco 
return 

A  1  =  1 .000 

RETURN 

ENG 


SUBROUTINE  B1V(M,P1) 
IMPLICIT  REAL *3 ( A -H  » 0- Z ) 
COMMON  A1 »E 1 ,A2  ,52 »° tINT 
GO  TO(l»2fi»2»2»l»2»l)*INT 
31=0SQRT(M*PI  ) 

R  ET'JR  N 

B1=0SQRT(  ( 2*M- 1  )*P1  ) 

RETURN 

ENO 
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SUBROUTINE  rT(X,Y,FX| 

IMPLICIT  RE  AL*3  (  a-'-i  ,Q-2  ) 

COMMON  a I ,R 1 , 42 ,52 ,R , I MT 
GO  T<?(  l  ,2 ,3,4,2 , 1 ,4,3  )  ,  I'lT 

F  X  =  Y*DS  I  N  (  2  .  ODQ*°':  X-X  )  /  GSO°  T  (  l.QOO-l.ODO/X) 
RETURN 

FX  =  Y'-DCGS  (2.0  00*°*X^X  ) /GSORT (  1 . OOJ- 1 . ODQ/X ) 
RETURN 

3  F  X  =  Y :X Q S I N  (  2 » 0 D 0 * 0 R- X R- X  )  / 0 5 G R  T  (  I  .000  *  l  ,000/X  ) 
RETURN 

FX=  YRQCOS  (  2.0D0*R*XRX  )  /  DSGRT  (  1  .COO*- 1  .000/X  ) 

RETURN 

END 


SU3RGUT I  ME  FO(X,Y,FX) 

IMPLICIT  REAL*? ( A-H,G-Z  ) 

COMMON  Al ,3L , A2,32  ,R,  IN T 
GO  T0(l,2, 3,4,1, 2*3,4), INT 

FX  =  CS  IN  (P.R-XMX*2  .GOO)  R-DSQRT  (  I.OOO/X-1 .000)  /  (  l  .OOO-Y/X  ) 
return 

F X  =0C OS  (  R*X*X*2 .000  ) *0$ DRT (  1 .000/X-l .000) / ( 1 .GOO-Y/X ) 
RETURN 

FX  =  DS IN (R*X*x*2.0DO ) *DSQRT(  I . 000/ X  - I . ODO ) / { I . 000  *  Y/X ) 
RETURN 

FX  =  OCO$ (R*X*X*2 .000) *QSQRT(  I.OOO/X-I .000) /( I .000* Y/X | 

RETURN 

END 
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